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ABSTRACT 

We present axisymmetric dynamical models of the edge-on SO galaxy NGC 4342. 
This small low-luminosity galaxy harbors, in addition to its outer disk, a bright nuclear 
stellar disk. A combination of observations from the ground and with the Hubble 
Space Telescope (HST) has shown that NGC 4342 rotates rapidly and has a strong 
central increase in velocity dispersion. 

We construct simple two-integral Jeans models as well as fully general, three- 
integral models. The latter are built using a modified version of Schwarzschild's 
orbit-superposition technique developed by Rix et al. and Cretton et al. These models 
allow us to reproduce the full line-of-sight velocity distributions, or 'velocity profiles' 
(VPs), which we parameterize by a Gauss- Hermite series. The modeling takes seeing 
convolution and pixel binning into account. 

The two-integral Jeans models suggest a black hole (BH) mass between 3 and 
6 X 10^ M0, depending on the data set used to constrain the model, but they fail 
to fit the details of the observed kinematics. The three-integral models can fit all 
ground-based and HST data simultaneously, but only when a central BH is included. 
Models without BH are ruled out to a confidence level better than 99.73 per cent. We 
determine a BH mass of S-Otl'l X Mq, where the errors are the formal 68.3 per 
cent confidence levels. This corresponds to 2.6 per cent of the total mass of the bulge, 
making NGC 4342 one of the galaxies with the highest BH mass to bulge mass ratio 
currently known. 

The models that best fit the data do not have a two-integral phase-space distribution 
function. They have rather complex dynamical structures: the velocity anisotropies 
are strong functions of radius reflecting the multi-component structure of this galaxy. 

When no central BH is included the best fit model tries to fit the high central 
velocity dispersion by placing stars on radial orbits. The high rotation velocities 
measured, however, restrict the amount of radial anisotropy such that the central 
velocity dispersion measured with the HST can only be fit when a massive BH is 
included in the models. 

Subject headings: black hole physics — galaxies: elliptical and lenticular, cD — 
galaxies: individual (NGC 4342) — galaxies: kinematics and dynamics — galaxies: 
nuclei — galaxies: structure. 
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1. Introduction 

Several lines of evidence suggest that active galactic nuclei (AGNs) are powered by accretion 
onto a super- massive black hole (BH) (Lynden-Bell 1969; Rees 1984). The much higher 
volume- number density of AGNs observed at redshift z ~ 2 than at z = 0, suggests that many 
quiescent (or 'normal') galaxies today must have gone through an active phase in the past, and 
therefore harbor a massive BH as well. Such a BH will significantly infiuence the dynamics of the 
galaxy inside a radius of influence, tbh = GMbh / , where cr is a characteristic velocity dispersion 
of the stars in the center. In particular, hydrostatic equilibrium requires that the rms velocities of 
the stars surrounding a massive BH follow an r~^/^ power-law (Bahcall & Wolf 1976; Young 1980). 

Since the late 70s, combined imaging and spectroscopy of the central regions of galaxies has 
suggested that massive BHs should be present in a number of early-type galaxies (see Kormendy 
&; Richstone 1995 for a review). Conclusive dynamical evidence for the presence of a central BH 
requires that a model with a BH can fit all observations (photometric and kinematic), and that no 
model without a BH can provide an equally good fit. Such conclusive evidence can only be inferred 
from observations that probe well inside the radius where the BH dominates the dynamics. Up to 
a few years ago, most claimed BH detections were based on observations with spatial resolutions 
of similar size as the radii of influence of the inferred BH masses (Rix 1993). This, together 
with the limited amount of freedom in the models used to interpret the data, has hampered an 
unambiguous proof for the presence of these BHs (i.e., the observed kinematics could not be 
confronted with all possible dynamical configurations without a BH). Often spherical models 
were used even when the observed fiattening was significant. If the models were axisymmetric, 
the distribution function (hereafter DF) was often assumed to depend only on the two classical 
integrals of motion, energy and vertical angular momentum; / = f{E,Lz). This implies that the 
velocity dispersions in the radial and vertical directions are equal (i.e., aji = az). It is well-known 
that strong radial anisotropy in the center of a galaxy results in a high central velocity dispersion, 
mimicking the presence of a massive BH (cf. Binney & Mamon 1982). Conclusive evidence for 
a BH therefore requires that one can rule out radial anisotropy as the cause of the high velocity 
dispersions measured, and models must thus be sufficiently general. 

Recently two major breakthroughs have initiated a new era in the search for massive BHs in 
normal galaxies. First of all, we can now obtain kinematics at much higher spatial resolution (down 
to FWHM ~ 0.1"), using specially-designed spectrographs, such as the Subarcsecond Imaging 
Spectrograph (SIS) on the Canada-France Hawaii Telescope, or the Faint Object Spectrograph 
(FOS) and STIS aboard the HST. This allows us to probe the gravitational potential much closer 
to the center, where the BH dominates the dynamics. Not only has this improved the evidence 
for massive BHs in several old BH-candidate galaxies (M31, Ford et al. 1998; M32, van dcr Marel 
et al. 1997, 1998; M87, Harms et al. 1994, Macchetto et al. 1997; NGC 3115, Kormendy et 
al. 1996a; NGC 4594, Kormendy et al. 1996b), but it has also provided new cases (M84, Bower 
et al. 1998; NGC 3377, Kormendy et al. 1998; NGC 3379, Gebhardt et al. 1998; NGC 4261, 
Ferrarese, Ford & Jaffe 1996; NGC 4486B, Kormendy et al. 1997; NGC 6251, Ferrarese, Ford & 
Jaffe 1998; and NGC 7052, van der Marel & van den Bosch 1998). Secondly, the revolutionary 
increase in computer power has made it possible to investigate a large number of fully general. 
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three-integral models based on the orbit-superposition method (Schwarzschild 1979). In the past 
deeade, this method has been used to build a variety of spherical, axisymmetric and triaxial 
models (e.g., Schwarzschild 1982; Pfenniger 1984; Richstone & Tremainc 1984, 1988; Zhao 1996). 
Levison &; Richstone (1985), Richstone & Tremaine (1985), and Pfenniger (1984) showed how to 
include rotation velocities and velocity dispersions as kinematic constraints. More recently, Rix et 
al. (1997) and Cretton et al. (1998) extended this modeling technique even further by fitting to 
the entire velocity profiles (see also Richstone 1997). Van der Marel et al. (1997, 1998) used this 
to build fully general, axisymmetric models of M32, and showed convincingly that M32 harbors 
a massive BH. Recent review papers on this rapidly evolving field include Ford et al. (1998), Ho 
(1998), Richstone (1998), and van der Marel (1998). 

In many galaxies where the presence of a BH has been suggested, a nuclear disk, seen close 
to edge-on, is present. These disks are either in gaseous form (M84, M87, NGC 4261, NGC 4594, 
NGC 6251, NGC 7052), or made up of stars (NGC 3115). It is easier to detect BHs in edge-on 
systems with disks, where one can use both the measured rotation velocities and the velocity 
dispersions to determine the central mass density. It is therefore not surprising that BHs have 
predominantly been found in galaxies with nuclear disks. Furthermore, nuclear disks allow a good 
determination of the central mass density of their host galaxies. Gaseous disks have the advantage 
that their kinematics can be easily measured from emission lines. Since gas in a steady-state 
disk can only move on non-intersecting orbits, the measured rotation velocities of a settled gas 
disk, in the equatorial plane of an axisymmetric potential, correspond to the circular velocities, 
Vc{R) = \/ i?d$ / di?. The rotation curve of a nuclear gas disk therefore provides a direct measure 
of the central potential gradient, and thus of the central mass density. However, often the gas 
disks are not in a steady state; many show a distorted morphology (e.g. M87, see Ford et al. 1994), 
and non-gravitational motion, such as outflow, inflow or turbulence can be present and complicate 
the dynamical analysis (e.g., NGC 4261, Jaffe et al. 1996; NGC 7052, van den Bosch & van der 
Marel 1995). Nuclear stellar disks do not suffer from this, but have the disadvantage that their 
kinematics are much harder to measure. First of all, the kinematics have to be determined from 
absorption lines rather than emission lines, and secondly, the line-of-sight velocity distributions, or 
velocity profiles (VPs), measured are 'contaminated' by light from the bulge component. However, 
van den Bosch & de Zeeuw (1996) showed that with sufficient spatial and spectral resolution 
one can resolve the VPs in a broad bulge- component and a narrow disk-component. From these 
VPs the rotation curve of the nuclear disk can be derived, providing an accurate measure for the 
central mass density. Therefore, galaxies with an embedded nuclear disk (either gaseous or stellar) 
observed close to edge-on are ideal systems to investigate the presence of massive BHs. 

In this paper we discuss the case of NGC 4342; a small, low-luminosity {Mb = —17.47) 
SO galaxy in the Virgo cluster. The galaxy is listed as IC 3256 in both the Second and Third 
Reference Catalogues of Bright Galaxies, since in the past it has occasionally been confused with 
NGC 4341 and NGC 4343 (see Zwicky & Herzog 1966). At a projected distance of ~ 30" SE of 
NGC 4342, a small galaxy is visible. It is uncertain whether this is a real companion of NGC 4342 
or whether it is merely close in projection. HST images of NGC 4342 revealed both an outer disk, 
as well as a very bright nuclear stellar disk inside ~ 1" (van den Bosch et al. 1994; Scorza & van 
den Bosch 1998). It is a normal galaxy, with no detected ISM (Roberts et al. 1991), and with small 
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color-gradients (van den Bosch, Jaffa & van der Marel 1998, hereafter BJM98). For its size and 
luminosity, it does however reveal a remarkably large central velocity dispersion and a very steep 
rotation-curve (see B JM98) . Unfortunately, the spectral resolution of the available kinematic data 
is insufficient to actually resolve the VPs in disk and bulge components. In order to determine 
the central mass density in NGC 4342, we thus have to construct dynamical models of the entire 
system: bulge and disk components. Here we present simple two-integral Jeans models as well as 
fully general three- integral models, and we provide evidence for the presence of a central massive 
dark object (MDO) of ~ 3 x 10^ M©. Throughout this paper we assume the MDO to be a BH, 
but we discuss alternatives in Section 7.2 

In Section 2 we briefly discuss the data used to constrain the models and in Section 3 we 
describe our mass model. In Section 4 we show the results of some simple two-integral modeling, 
and we discuss its shortcomings. Section 5 describes the general outline of the three- integral 
modeling technique. In Section 6 we discuss shortcomings of the velocity-profile parameterization 
used when applied to dynamically cold systems, and present a modified approach. The results 
of the three-integral modeling are discussed in Section 7. Finally, in Section 8, we sum up and 
present our conclusions. Throughout this paper we adopt a distance of 15 Mpc for NGC 4342, 
consistent with the distance of the Virgo cluster (Jacoby, Ciardullo & Ford 1990). 



2. The data 

All data used in this paper are presented and discussed in detail in BJM98. Here we merely 
summarize. 



2.1. Photometric data 

BJM98 used the Wide Field and Planetary Camera 2 (WFPC2) aboard the HST to obtain 
U, V and / band photometry of NGC 4342. The spatial resolution of these images is limited by 
the HST Point Spread Function (FWHM ~ 0.1") and the size of the pixels (0.0455" x 0.0455"). 
The full field-of-view covers about 35" x 35". Figure || shows contour plots at two different scales 
of the I-band image. The presence of the nuclear disk is evident from the highly flattened, disky 
isophotes inside 1.0". 



2.2. Kinematic data 

Using the ISIS spectrograph mounted at the 4.2m William Herschel Telescope (WHT) at La 
Palma, BJM98 obtained long-slit spectra of NGC 4342 along both the major and the minor axis. 
The spectra have a resolution of cxinstr = 9kms~^, and were obtained with a slit width of 1.0" under 
good seeing conditions with a PSF FWHM of 0.80" (major axis) and 0.95" (minor axis). After 
standard reduction, the parameters (7, y, a, /i3, /14) that best fit the VPs were determined using 
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the method described in van der Marel (1994). These parameters quantify the Gauss-Hermite 
(GH) expansion of the velocity profile C{v) as introduced by van der Marel & Franx (1993): 

7 



where 



and 



C{v) = j-a{w)\l + Y,h,Hj{w)\, (1) 

i=3 



w^{v- V)/a, (2) 



a[w) = ^e-i-^ (3) 



Here v is the line-of-sight velocity, Hj are the Hermite polynomials of degree j, and hj are 
the Gauss-Hermite coefficients. The first term in equation ([l|) represents a Gaussian with line 
strength 7, mean radial velocity V , and velocity dispersion a. The even GH-coefficients quantify 
symmetric deviations of the VP from the best-fitting Gaussian, and the odd coefficients quantify 
the asymmetric deviations. 

We have averaged the kinematic WHT data at positive and negative radii. In this way we 
obtain sets of (y, o", /13, /i4) at 19 different positions along the major axis and 8 along the minor 
axis. 

BJM98 also obtained FOS spectra at 7 different aperture-positions, all inside the central 0.5" 
of NGC 4342, using the circular 0.26"-diameter aperture (the FOS . 3-aperture). Due to the 
limited signal-to-noise ratio (S/N) of these spectra, only (7, V, a) of the best-fitting Gaussian could 
be determined. Rotation velocities and velocity dispersions along the major axis of NGC 4342 
for both the WHT (crosses) and FOS (solid dots) are shown in Figure The central rotation 
gradient, as measured with the FOS, is extremely steep {V ~ 200 km s^^ at 0.25" from the center). 
In addition, the velocity dispersion increases from ~ 90 kms~^ at the outside (the 'cold' outer 
disk) to 317kms~^ in the center as measured with the WHT. The central velocity dispersion 
increases to 418 kms~^, when observed at four times higher spatial resolution with the FOS. 



3. The mass model 

We have used the Multi-Gaussian Expansion (MGE) method developed by Emsellem, Monnet 
& Bacon (1994, hereafter EMB94) to build a mass model for NGC 4342. The method assumes 
that both the PSF and the intrinsic surface brightness are described by a sum of Gaussians, 
each of which has 6 free parameters: the center {xj,yj), the position angle, the flattening q'j, the 
central intensity /j, and the size of the Gaussian along the major axis, expressed by its standard 
deviation a'j. The best-fitting parameters of the different Gaussians are determined using an 
iterative approach in which additional components are added until convergence is achieved (see 
EMB94 for details on the method). This method is well suited for complicated, multi-component 
galaxies such as NGC 4342. 
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We fitted the HST /-band PSF by a sum of 5 circular (i.e., q'^ = 1) Gaussians (see BJM98). 
Using this model PSF we derived the parameters of the N Gaussians describing the intrinsic 
surface brightness (i.e., deconvolved for PSF effects) by fitting to the HST /-band image of 
NGC 4342. We forced the A'^ Gaussians to have the same position angle and center, which yields 
an axisymmetric mass model (see below). Therefore, the model is described by 3A^ + 3 free 
parameters, which are simultaneously fit to the image. We achieved convergence with = 11 
Gaussian components. The results of the fit are shown in Figure |^, where we show contour plots 
of the /-band image with superimposed contours of the convolved surface brightness of the MGE 
model. The fit is excellent, except for a small discrepancy at the outside. This is due to slight 
twisting of the isophotes at large radii (see BJM98). Since our model is axisymmetric , this cannot 
be modeled. Nevertheless, the discrepancy is small, and is unlikely to affect our conclusions on the 
dynamics of the central region. The parameters of the different Gaussian components are listed in 
Table 1. 



The total luminosity of the MGE model in the /-band is L/ = 3.57 x 10^ Lq. This 
yields M/ = —19.86. The absolute blue magnitude of NGC 4342, at a distance of 15 Mpc, is 
Mb = —17.47 (Sandage & Tammann 1981), and we thus find B — I = 2.39. This is consistent 
with the colors of NGC 4342 presented by BJM98. They find J7 - F « 1.5 and F - / « 1.3. 
We thus derive B — V k, 1.09, in good agreement with the average value for early-type galaxies 
(Faber et al. 1989). If we assume that the luminosity distribution of the bulge corresponds to 
the Gaussian components rounder than q'- = 0.3, we find that the bulge makes up ~ 52 per 
cent of the total luminosity of NGC 4342. The outer disk, described by Gaussian components 9 
and 10, makes up an additional 46.5 per cent, and the nuclear disk (modeled by Gaussian 
component 4) adds only about 1.5 per cent to the total luminosity. There is no reason that the 
mathematical components correspond to actual physical components, but at least this gives an 
order-of-magnitude description of the luminosities of the bulge and the two disk components. A 
more accurate disk-bulge decomposition, which yields similar results, is discussed in Scorza & van 
den Bosch (1998). 

Assuming that the density is built up from a sum of three-dimensional Gaussians stratified on 
spheroids, one can, for any inclination angle i, analytically calculate the density distribution from 
the MGE fit to the intrinsic surface brightness. The mass density of such an MGE model is given 
by 



p{R, z) = T ^ /j-exp 



1 



+ 72 



(4) 



where T is the mass-to- light ratio, and Ij, aj and qi are related to /j, a'j, q'j and i (see EMB94). 
The potential that corresponds to this density distribution follows from solving the Poisson 
equation. This yields 



$(/?, z) = -inGT a^Qjlj exp 



2aj 



b;' 



dt 



(5) 



where = I — q^. 
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The inclination angle is well constrained by the thinness of the nuclear disk: i > 83deg 
(Scorza & van den Bosch 1998). Throughout we assume that NGC 4342 is observed edge-on (i.e., 
i = 90deg). Given the lower limit on the inclination angle of 83deg, this assumption does not 
significantly influence the conclusions presented in this paper. 



4. Jeans modeling 
4.1. Formalism 

The three-integral modeling described in the next section requires large amounts of CPU 
time. We therefore decided to first explore parameter space of the models (i.e., mass-to-light ratio 
and mass of the possible BH) by solving the Jeans equations and assuming that the phase-space 
distribution function depends only on the two classical integrals of motion. The Jeans equations 
for hydrostatic equilibrium are moment equations of the collisionless Boltzmann equation (see 
Binney & Tremaine 1987). They relate the velocity dispersion tensor and the streaming motion 
V to the density p and potential For an axisymmetric system with distribution function 
f{E, Lz) one always has ur = az and vrVz = 0, and the Jeans equations in cylindrical coordinates 
reduce to 

*«'+pf^M + l)=o, (0) 



dR R dR 

^+pf=0, (7) 

Here 7 denotes the local average over velocities. From equation (0) and p(R, z) and ^{R,z), one 
can, at every point {R, z) in the meridional plane, calculate a\ (= a1) by simple integration. By 
rewriting equation (P), one can compute = v^'^ + cr^ without the need of performing a numerical 
(ill-conditioned) derivative (see e.g. Hunter 1977; Simien, Pellet & Monnet 1979; Binney, Davies 
& Ilhngworth 1990). 

The expressions for a\ and for the density distribution of equation ^ are given in EMB94 

(their equations [42] and [44]). The Jeans equations do not prescribe how splits in streaming 
motion and azimuthal velocity dispersion cr,^. We follow the approach introduced by Satoh 
(1980), and write 



'^^ = k\/vl-a% (8) 

such that we can control the anisotropy cr(j>/crfi by means of the free parameter k. For k = 1 
the model is fully isotropic with a,^ = ctr = cr^. Once k has been fixed one can project the 
luminosity- weighted dynamical quantities on the plane of the sky {x',y') to yield the projected 
rotation velocities 

— 1 f°° 

yvot{x',y') = — -r i^v^smicos(j) dz', (9) 

and the rms velocities 

K^msix', y') = -Wf-; — j\ I ^{^z '^o^^ * + '^R (^sin^ i + (u^ + U^^) cos^ 0sin^ i) dz'. (10) 
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Here u = p/T is the luminosity density and S{x',y') is the projected surface brightness at position 
{x',y'). Vrot and Vrms are the true first and second order moments of the hne-of-sight velocity 
distribution. The projected velocity dispersion ap{x',y') is simply derived from 

^pix', y') = ^Jv,l,,{x',y')-VUx',y'). (11) 

It is straightforward to include a BH in such a model, by simply adding GMbh/V-R^ + z"^ to the 
stellar potential (^). 

4.2. Application to NGC 4342 

We use the Jeans equations to calculate the predicted rotation velocities and velocity 
dispersions for the luminosity distribution of NGC 4342. We assume that the stellar mass-to-light 
ratio T and the anisotropy parameter k are constant throughout the galaxy. We calculate V-^ot and 
K-ms (using equations Q and |jl^) on a two-dimensional grid on the sky. The grid is logarithmic 
in r (in order to properly sample the strong gradients near the center), and linearly sampled in 0. 
Once Vrot and T^^g are tabulated, we convolve them with the PSF of the observations, weighted 
by the surface brightness. After pixel binning, taking the proper slit width into account, these are 
compared to the observations. 

The V and a determined from the GH fitting to the WHT VPs cannot be compared directly 
to V^ot and dp derived from the modeling discussed above: the latter ones correspond to the 
true moments of first and second order of the VPs, whereas the former ones correspond to the 
best-fitting Gaussian. We therefore recalculated the VP from V ^ a, /13 and /i4, from which we then 
estimate the first and second order moments for direct comparison with the Jeans models. 

The results are shown in Figure |3|, where we plot V^ot, K-ms and ap of the VPs along the major 
axis for both the WHT and the FOS data. Also plotted are predictions for four models, that 
only differ in the mass of the central BH (0, 3, 5 and 10 x 10^ M©). All models have i = 90deg, 
T/ = 6.2 M0/ Lq and k = 1 (i.e., all models are fully isotropic). For this value of k, we obtain the 
best fit to the observed velocity dispersion outside ~ 2". However, 14ot is not very well fitted: the 
wiggles in the rotation curve are not reproduced by the model. One can alter k, as function of 
radius, such that we fit these wiggles, but at the cost of introducing them in the velocity dispersion 
profile. This is due to the poor fit of the Jeans models to the rms velocities (see lower-left panel 
of Figure ^). depends only on the sum of cj^ and v^"^ (see equation [jl^) and is therefore 

independent of k. Consequently, the Jeans models cannot simultaneously fit T^ot and ap along 
the major axis. This suggests that the assumption made, i.e., / = f{E,Lz), is wrong, and that 
three-integral models are required. 

The Jeans models without a central BH clearly underpredict the central velocity dispersion 
(for both the WHT and the FOS measurements), as well as the central rotation gradient measured 
with the FOS. Models with a massive BH provide a much better fit. The actual mass of the BH 
depends on the data set used to constrain the model: the WHT data suggests a BH mass of 
~ 3 X 10^ Mq, whereas the FOS data are best fitted with Mbh « 6 x 10^ Mq. So although the 
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two-integral Jeans modeling cannot fit all the observed kinematics, it does suggest that a BH of a 
few times 10^ M0 may be present in the center of NGC 4342. 

In Figure ^ we plot the ratio of the observed rotation velocity over the rotation velocity of 
the best-fitting isotropic Jeans model, V^bs/V"mod (along the major axis), versus the local observed 
ellipticity (in /-band) of NGC 4342. There is a clear correlation in the sense that the Jeans 
models underpredict the rotation velocity in the strongly flattened region, and overpredict Vrot 
in the less flattened region. The ellipticity is a measure of the local disk-to-bulge ratio, and 
this therefore suggests that disk and bulge have different velocity anisotropics. Since Vrot scales 
with \/T, another possibility may be that the disk and bulge are made up of different stellar 
populations (whose mass-to-light ratios are different by almost a factor two). However, the 
separate components (bulge, nuclear disk and outer disk) do not stand out as separate entities in 
either the U — V 01 the V — I color images (see BJM98), rendering this explanation improbable. 



5. Three-integral modeling 

In order to investigate the presence of a ~ 3-6 x 10^ Mq BH as suggested by the Jeans 
models, we now construct fully general, axisymmetric models of NGC 4342. We use an extension 
of Schwarzschild's (1979) orbit-superposition method (see de Zeeuw 1997), developed by Rix et 
al. (1997, hereafter R97) and Cretton et al. (1998, hereafter C98). The main method, but for 
spherical systems, is outlined in R97. The application to axisymmetric systems is discussed in 
C98. Here we briefly outline the method, and we refer the interested reader to R97 and C98 for 
details and tests of this modeling technique. 



5.1. The method 

The first step of the method is to integrate orbits in the combined potential <I>stars + ^bh- 
Each orbit is then projected onto the space of observables, taking convolution with the PSF and 
pixel binning into account. Finally, a non- negative least-squares algorithm is used to determine 
the distribution of orbit weights that best fits the observational data (taking the observational 
errors into account), while also reproducing the luminous density distribution of the model. 

Throughout we limit ourselves to models with an inclination angle i = 90deg, and we 
assume that the stellar population has a constant mass-to-light ratio. Therefore, each model is 
characterized by only two free parameters: the mass-to-light ratio, T/, and the mass of the black 
hole, Mbh- Our aim is to find the set (Mbhj '^/) that best fits the available constraints (surface 
brightness and velocity profiles). 
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5.2. The orbit library 



The motion of a star in an axisymmetric potential, for which E and = Rv^ are 
conserved, can be reduced to motion in the meridional (i?, 2)-plane, in the effective potential 
$e = ^{R, z) + L1/{2R?). The orbit is constrained within a region bounded by the zero velocity 
curve (ZVC) defined through E = ^^{R,z). 

Each orbit in an axisymmetric potential z) admits two integrals of motion: energy 

E = ^{R, z) + ^ and vertical angular momentum = Rv^. Regular orbits admit one additional 
integral, I3, which in general is not known analytically. Such an orbit is confined to a sub-space 
inside the ZVC. We only found a very small fraction of our orbit catalog to be irregular. 

The orbit library has to be set up such that one properly samples the full extent of phase 
space. The sampling has to be sufficiently dense to suppress discreteness noise, but is limited by 
the amount of available CPU time. After some testing we have chosen to calculate orbits on a 
20 X 20 X 7 {E, Lz, l3)-gnd. Each energy is uniquely defined by a circular radius Rc according to 



The 20 energies are sampled logarithmically between Rc{E) = 0.01" and Rc{E) = 60". These 
values were chosen such as to encompass the major fraction of the total mass of the galaxy. The 
mass inside Rc = 0.01" is only a fraction of 3.25 x 10~^ of the total mass, whereas a fraction of 
~ 2.85 X 10~^ of the mass is located outside Rc = 60". For each energy, we calculate the maximum 
vertical angular momentum 



which corresponds to the circular orbit with energy E^ and sample Ir/I = \Lz\/ Lz^max on a linear 
grid of 10 values between 0.01 and 0.99. Hence, the purely circular and radial orbits are presumed 
to be represented by their closest neighbors on the grid, but are not explicitly calculated. At 
each value of I77I, only the orbit with positive angular momentum has to be integrated, since its 
counterpart (L^ = —rjLz^max) is simply a mirror refiection around zero velocity. Since the third 
integral I3 can generally not be expressed explicitly in terms of the phase-space coordinates, we 
use the method suggested by Levison & Richstone (1985), and take the starting point on the 
ZVC as a numerical representation of the third integral (see also C98). For that purpose we first 
calculate, for each (E, Lz)-paiT, the locus (i?tt>-Ztt) on the ZVC of the 'thin tube orbit'. For a 
certain value of E and L^, this is the only orbit that touches the ZVC at only one value of R. 
All other regular orbits touch the ZVC at two different values of R. We sample I3 by linearly 
sampling the ii-coordinate of the orbit's starting point on the ZVC (Rzvc) between R^ and -Rmax, 
where -Rmax is the maximum extent of the ZVC in the equatorial plane. 

Since the calculation of the potential (equation ) and the forces all require the evaluation 
of a numerical quadrature, we have calculated them on a 4000 x 300 {R, 9)-grid in the meridional 
plane. The grid is sampled linearly in 9 between and vr/2, and logarithmically in R between 
10~^ and 10'^ arcseconds. Each orbit is integrated for 200 radial periods, using linear interpolation 




R—Rc 



+ ^(i?c,0). 



(12) 




(13) 
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between grid points to evaluate the potential and forces. On average, the energy conservation 
over 200 radial periods is better than one part in 10^, justifying the interpolation scheme adopted. 
In total 20 X 10 X 7 = 1400 orbits are integrated, resulting in a library of 2 x 1400 = 2800 orbits 
(when doubling for the — orbits). 

During the integration of each orbit we project its phase-space coordinates onto the space 
of observables (x', y', vios), where {x',y') is the plane on the sky, and fios is the line-of-sight 
velocity. In the following we use v as shorthand for v\os- We adopt a three-dimensional grid in the 
(x', y', w)-space, i.e., our storage cube, in which we record the fractional time the orbit spends in 
each of the cells (see R97 for details). Once the orbit integration is finished, we convolve each of 
the velocity slices (x', y') with the PSF appropriate for the observations. Since we have kinematic 
constraints obtained with three different instrumental setups and different PSFs (WHT major 
axis, WHT minor axis, and FOS apertures) we use three separate (x', y', f )-storage cubes each 
of which is convolved with the respective PSF. We use two cubes with 0.1" x 0.1" (x', y')-cells 
for the WHT major and minor axes. For the FOS-cube, 0.05" x 0.05" cells are used to comply 
with the higher spatial resolution of the HST. For all storage cubes we use 101 velocity bins of 
30 kms~^. The final step is to calculate the contribution of the orbit to each of the positions in 
the plane of the sky where we have photometric and/or kinematic constraints. These positions are 
in general extended areas (e.g., determined by the pixel size of the CCD, the slit width and the 
pixel rebinning used to obtain spectra of sufficient S/N). For each constraint position / = 1, Nc 
we therefore sum the fractional times over the area of / (see R97 for details). This gives us in the 
end, for each orbit k and each constraint position /, the properly PSF-convolved velocity profile 
'histogram' VPf^, integrated over the area of position /. By using fractional times we ensure that 
each orbit is normalized to unity. 



5.3. The observational constraints 

The final step of the orbit-superposition method is to find the set of non-negative weights, 
7fc, of each orbit that best matches the kinematic constraints and reproduces the luminous mass 
density of the model. Since we normalize each orbit to unity, the orbital weights 7^ measure 
the fraction of the total light of the galaxy that resides on orbit k. We use the following sets of 
constraints: 

The solution 7 has to reproduce the luminosity density p{R,z)/T (equation [^). We have 
subdivided the first quadrant of the meridional plane in 20 x 7 (r, 0)-cells, with r and 9 the 
standard spherical coordinates. The grid that encompasses those cells is binned linearly in 9, and 
logarithmically in r between 0.01" and 60". For each orbit k we store the fractional time spent 
in cell n. Since we integrate the orbits in the meridional plane, this is similar to the fractional 
time the orbit spends in the three-dimensional volume obtained by integrating the area of the 
cell over 27r radians in the (/>-direction. Therefore we have computed the total luminosity, of 
the MCE model inside each volume n. We will refer to these constraints as the "self-consistency 
constraints" . 
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From the WHT spectra we obtained sets of {Vi, ai, h^^i, h^^i) at 19 positions I along the major 
axis, and 8 along the minor axis. The quantities are the luminosity- weighted averages over the 
areas of constraint positions I = 1,...,27. For each of these 27 positions we have calculated the 
surface brightness Si integrated over that area, and convolved with the appropriate PSF. The 
quantities Si, V/, (T;, h^^i and h^^i parameterize the velocity profiles £f^^{v) with 

/oo 
Cf\v)dv. (14) 

One can rewrite this parameterization in the form (5^, Sihi^i, Sih2,i, Sih^^i, Sih^^^i), with 

/oo 
Cf\v)a{wi)Hm{wi)dv, (15) 
-oo 

where m = 1, 4, a is again the standard Gaussian (see equation Q), and 

wi = , (16) 

with Vi and ai the measured rotation velocity and velocity dispersion of the VP's best-fitting 
Gaussian at constraint position /. Note that with this definition, hi^i and /i2,z are zero for 
all constraint positions. This parameterization has the advantage that the orbit-superposition 
problem for the orbit weights is linear (see R97). 

In addition, we obtained (V;,cr;) at 7 FOS aperture-positions. Again we have calculated, for 
each of these positions I, the PSF-convolved, aperture-integrated surface brightness Si. The PSF 
used and the aperture diameter adopted are described in BJM98. As for the WHT spectra, we use 
the parameterization [Si, Sih^^i with m = 1,2) as constraints rather than (S;, V/, cr/). 

In total we thus have 140 self-consistency constraints 34 constraints on projected surface 
brightness S^ and 122 kinematic constraints Sih^i. 



5.4. Non- negative least squares fitting 

At each constraint position I we have a measured velocity profile Cf^^{v), which we 
parameterized by {Si, Sihm.i)- For each of the observational constraints, we have the measurement 
errors ASi, AVi, Aai, Ah^^i and Ah/^^i from which we can calculate the errors A{Sihmj)- As 
described in Section 5.2 we determined, for each orbit k and each position I, the velocity profile 
VPf^ at velocity v. We parameterize each orbital VP by {Sf,Sfh^i) with m = 1,...,4 using 
equations (14) and ([Tsl ) and with £f^^{v) replaced by VPf^, and by changing the integration to a 



summation over all velocity bins. The orbit weights {k = 1, A^orbits) that result in the best 
fit to the observations can be determined by minimizing 




^{Slhml) 



(17) 
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In addition to minimizing Xobs' want the solution to match the luminosity density in the 

meridional plane, i.e., we also want to minimize 



where AL„ sets the accuracy for reproducing the luminosity density in the meridional plane. 
Throughout, we set ASi = 0.015; and ALi = O.OIL;, such that we aim for an accuracy of one 
per cent in reproducing both the projected PSF-convolved surface brightness and the luminosity 
density in the meridional plane. It is in principle sufficient to fit only the luminous density in 
the meridional plane: the surface brightness should be fitted automatically. However, because of 
discretization this is in practice not necessarily so, and we thus include the surface brightness at 
the constraint positions I as separate constraints. 

Minimizing = Xobs + Xsc * least-squares problem for a matrix equation (see R97 for a 
detailed description of the matrices involved). It has to be solved under the physical constraint 
that 7fc > 0. Following Pfcnniger (1984), R97 and C98, we use the Non-Negative Least Squares 
(NNLS) algorithm by Lawson &i Hanson (1974) to solve for the orbit weights. 



The NNLS matrix equation solved is numerically rather 'ill-conditioned' giving rise to a 
distribution function with strong oscillatory behavior. Such DFs are unphysical (e.g., Lynden-Bell 
1967; Spergel & Hernquist 1992; Merritt 1993). Smoothing in the solution space can be achieved 
via rcgularization (e.g., Merritt 1996 and references therein). We follow the scheme used by Zhao 
(1996), which is based on a minimization, up to a certain degree, of the differences in weights 
between neighboring orbits. The technique is described in R97, C98, and van der Marel et al. 1998, 
and we refer the interested reader to those papers for details. The extra rcgularization constraints 
result in a less good fit to the data. The amount of smoothing is set such that the regularized 
model is still compatible with the data in a statistical sense, i.e., such that Ax^ = ~ Xmin ~ 1' 
with Xmin the value obtained without rcgularization. Unless mentioned otherwise, we discuss 
results in which no rcgularization has been adopted. 



Although the modeling technique outlined in the previous sections works well for dynamically 
hot galaxies (e.g. M32, see van der Marcl ct al. 1998), some problems arise when trying to apply it 
to dynamically cold systems such as NGC 4342. This is due to the fact that we do not include the 
zeroth-order moment ho in the fit. This quantity measures the normalization of the best-fitting 
Gaussian to the normalized VP and is observationally inaccessible: it is directly proportional to 
the unknown difference in line strength between the galaxy spectrum and the template spectrum 
used to analyze. (In practice one uses the assumption /iq = 1 to estimate the line strength from 




(18) 



5.5. Regularization 



6. Shortcoming of the Gauss-Hermite parameterization 
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the observations.) As we now show, excluding /iq from the GH series and expanding the orbital 
VPs around the observed VPs up to fourth order only, can lead to artificial counter-rotation in 
models of cold systems. 

Assume we have an observed velocity profile VPobsi at a certain position on the sky, that is 
perfectly Gaussian. The GH moments hm with m > of such a VPobs will all be zero. In the 
method described above, we derive the orbital /i,^ (m = 1, ...,4) from equation ( p^ ) in which the 
observed V and a enter in the weighting function a{w)Hm{w). In the NNLS algorithm we solve 
for the orbit weights by minimizing the difference between the GH moments of VPobs and VPorb- 
In principle, the differences between the Sih^.i and YliklkSfh^i are minimized, but for simplicity 
we illustrate the problem with a one-orbit model. 

In the ideal case, an orbit whose VP deviates more strongly from the observed VP will be 
assigned a smaller weight. However, this is not always the case with the VP parameterization 
described in Section 5.3. To illustrate this we calculated the GH moments {m = 0, ...,2) of 
a Gaussian expanded around another Gaussian, both with the same dispersion a. In Figure |5| 
we plot the resulting moments as function of the velocity difference AF between the two 
Gaussians (in units of a). For Al/ = 0, /iq = 1 and all higher-order moments are zero (i.e., the 
two Gaussians are identical). In the regime |Ay|/cj <^ 2 the higher-order moments increase with 
increasing velocity-difference. For larger values of |Ay|/cj they start to decrease again to reach 
approximately zero for |Ay|/cj ^ 5. On the contrary, /iq equals unity when the VPs are identical, 
and decreases monotonically for increasing |Ay|/o". 

In NGC 4342, which is a dynamically cold system, the major-axis kinematics reach 
V^bs/cobs > 2.5 at the outside (see Figure |6|). An orbital VP with approximately the same V and a 
as an observed VP with V/a > 2.5 will have its hi and /12 close to zero. Consequently, it will likely 
be given a non-zero weight. The same orbit, but with opposite sense of rotation (i.e., with reversed 
vertical angular momentum), will also have hi and /12 close to zero, since |Ay|/cj > 5. In other 
words, the +Lz and —L^ orbits have the same S, hi and /12 and are therefore indistinguishable 
for the fitting algorithm. The +Lz and —L^ orbits have opposite values of /13, but if |/i3| is small, 
the difference between the two orbits, in terms of the VP parameterization used, remains small. 
If it would have been possible to include /iq in the fit, then the +Lz would have /iq close to 
unity, and the —L^ orbit would have /iq ~ (see Figure ^: the two VPs would have been easily 
distinguishable. In addition, if the observational data had been of sufficient S/N to derive GH 
moments up to very high order, it would also have been possible to distinguish the VPs of the 
+Lz and — orbits. In practice, however, one requires unrealistic high S/N spectra to be able to 
measure these moments. 

We can thus expect that our solutions will have significant amounts of counter-rotation at 
radii where V^bs/^'obs ^ 2. We indeed found solutions in which the reconstructed VPs at large 
radii along the major axis have two peaks at positive and negative vios (see Section 7.1). Van 
der Marel et al. (1998) applied the same modeling technique to the dynamically hot system M32, 
which has V/a < 1.0 everywhere, and therefore did not encounter this problem. 

To solve the problem outlined above, we use a 'modified approach' in which we add some extra 
constraints to the models: we exclude counter-rotating orbits whose circular radius Rc is larger 
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than a limiting radius Rum- We have chosen Rn^ = 4.5", since for larger radii K3bs/o"obs > 1-5 (see 
Figure and the VP parameterization used becomes poor in the light of the problem discussed 
above. As we show in Section 7.1, this solves the problem of the artificial counter-rotation at 
large radii along the major axis, but it has the disadvantage that our models are not fully general 
any more: we have imposed some a priori constraints on the amount of counter-rotating orbits 
in NGC 4342. On the other hand, the actual observations do not reveal any counter rotation at 
large radii. Critically, one may argue that the VP analysis of the spectra is not suitable to detect 
such counter rotation, since we do not go to sufficient high order in the GH expansion. We have 
therefore also analyzed our spectra with the unresolved Gaussian decomposition (UGD) method 
(Kuijken & Merrifield 1993), which is not hampered by this limitation (e.g., Merrifield & Kuijken 
1994), and found good agreement with the results of the GH parameterization. This suggests 
that indeed NGC 4342 has a negligible amount of counter-rotating stars in its outer disk. The 
additional constraint imposed on the models, when using the modified approach, is thus justified 
observationally. 



7. Results and discussion 

7.1. The black hole and mass-to-light ratio 

Based on the BH mass and mass-to-light ratio suggested by the Jeans modeling (Section 4.1) 
we construct three-integral models with M^u/Tj in the range — 1.5 x 10^ Mq and 4 < T/ < 9. 
For each value of M^n/Tj only one orbit-library has to be constructed: a change in mass-to-light 
ratio is equivalent to a scaling of the model velocities proportional to -v/T- We therefore calculate 
10 different orbit libraries, all with Tj = 1, that differ only in the mass of the central BH. If one 
scales to another value of T/, the mass of the BH changes accordingly to T/ Mbh- We sample the 
mass-to-light ratio at 16 different values in the interval T/ G [4,9]. 

For each of the in total 160 different (Mbh, T"7)-niodels we minimize = Xobs + Xsc 
(equations |jl^ and |jl8|). We use statistics to compare different (Mbh, Tf'/)^models in a proper 
statistical way. We determine the measure Ax^ = ~ Xmin' where Xmin overall lowest x^- 

Under the assumptions that the errors are normally distributed and that there are no numerical 
errors in the model, one can assign confidence levels to the measure A^^. The exact level of 
confidence depends on the number of degrees of freedom in the models (Press et al. 1992), in our 
case two: Mbh and T/. 

The resulting x^-plots are shown in Figure 0. The first three contours show the formal 68.3, 
95.4 and 99.73 per cent confidence levels (the latter one is plotted with a thick contour). The solid 
dots correspond to actual model calculations. Bi-cubic spline interpolation is used to calculate 
at intermediate points. Four different plots are shown, labeled a to d. Plot a (upper-left panel) 
shows the plot of the fits, when using only the V and a measurements of the WHT spectra as 
constraints. Clearly, a large range of parameter space gives equally good fits. There is at best only 
a marginal indication that models with a BH fit the ground-based rotation velocities and velocity 
dispersions better than without a BH. Plot b (upper-right panel) is similar to plot a, except that 
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we have now included the /13 and /14 measurements of the WHT spectra as constramts on the 
models. Two changes are evident. First, the allowed range of T/ at given BH mass is much smaller. 
Secondly, BHs more massive than ~ 6 X 10^ Mq can be ruled out at the 99.73 per cent confidence 
level. Plot c (lower-left panel) shows the resulting x^-plot when all constraints, including the HST 
measurements are taken into account. These high spatial-resolution measurements allow us to 
rule out models without a BH at a confidence level better than 99.73 per cent. Finally, in plot 
d (lower-right panel) we show the results when using the modified approach on the entire set of 
constraints. As can be seen, the exclusion of counter-rotating orbits with circular radii beyond 
4.5" puts some further limits on the allowed range of acceptable models, consequently contracting 
slightly the 99.73 per cent confidence region. 

Although the solutions we find with the standard approach can result in models with 
significant counter-rotation (see Section 6), one can still meaningfully use the statistics to put 
confidence levels on Mbh and T/. The main requirement is that orbital VPs and observed VPs 
are parameterized in exactly the same way, which is the case. The counter rotation that we find 
is simply a consequence of our particular VP parameterization. Although the modified approach 
results in a more strict solution-space, and in principle is based on additional constraints that 
are observationally justified, we will nevertheless consider the surface of Figure |^ as the main 
result. We merely present the results of the modified approach to show that the counter-rotation 
in the best-fitting models has no significant infiuence on the BH mass and mass-to-light ratio of 
our best-fitting model. Furthermore, Figure |^ results in more conservative estimates of the errors 
on BH mass and mass-to-light ratio. 

The labeled asterisks in panel c indicate three models that will be discussed in more detail 
below. Model B provides the best overall fit and has Mbh = 3.6 x 10® Mq and T/ = 6.25 Mq/ Lq. 
From the statistics we derive Mbh = 3.0+{;o x 10® Mq and T/ = 6.3^0:4 Mq/ Lq. The errors 
are the formal 68.3 per cent confidence levels. For comparison, the modified approach yields 
Mbh = S.O^ol x 10® M© and T/ = 6.5to| Mq/ Lq, which is consistent with and shghtly more 
strict than the results obtained with the standard method. 

Van der Marel et al. (1998) used the same three- integral modeling technique and statistical 
means to infer that M32 harbors a massive BH of (3.4 it 0.7) x 10^ Mq. To facilitate direct 
comparison, the surfaces presented here are plotted with the same contour levels as in van der 
Marel et al. (1997, 1998). The contours of Figure correspond to the goodness-of-fit to all 
constraints. Van der Marel et al. plotted contours of x^ of the fit to the kinematic constraints only 
(i.e., the right term in equation 0]), but noted that contour plots of the total x^ look similar. 
We have computed the surfaces of the x^ that only measures the fit to the kinematic constraints 
and indeed found very similar contours: all models can fit the surface brightness and luminosity 
density to better than one per cent as requested. 

In order to compare the kinematical predictions of the models with the actual observations, 
we compute the velocity profile of the model at each constraint position /, VPi^v, as the weighted 
sum of all VPf^, using the solution 7 for the orbit weights. For these VPs we then compute Vi 
and ai of the best-fitting Gaussian, as well as 5/ and the GH moments h^^i and h^^i, all of which 
can be directly compared to the observations. The kinematical predictions of models A, B and C 
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are plotted in Figure ^ together with the observations. Although we did not measure the /13 and 
/i4 coefficients for the FOS spectra, we plot the model predictions for these quantities. All three 
models provide equally good fits to the projected surface brightness and the meridional luminosity 
density (not plotted). The main difference between models A and B is their fit to the central HST 
velocity dispersion. The model without BH (model A) underpredicts the observed dispersion by 
~ 130 kms~^ (~ 4.OA0"). The rapid central rotation in NGC 4342, as measured with the FOS 
is remarkably well fit even without central BH. The main difference between models B and C 
is their fit to the central WHT velocity dispersion along the minor axis, and the fit to the /13 
measurements along the major axis. Although model C fits the central HST velocity dispersion 
even better than model B, the fit to the higher-order GH coefficients is worse to such an extent 
that this model can be ruled out to better than 99.73 per cent confidence (cf. panels a and b of 
Figure |^ . 

The large differences between the predictions of model C at the outermost point along the 
major axis (WHT measurements) and the actual observed values is due to the problem with 
the counter-rotation. This is illustrated in Figure where we plot VPs at four different radii 
along the major axis. Solid dots represent the 'observed' velocity profiles, reconstructed from the 
measured GH-parameters (assuming that all GH coefficients of order five and higher are zero). 
The solid lines represent the reconstructed VPs of model C, and the dashed lines correspond to 
the Gauss-Hermite series fitted to these model VPs. The upper four panels show the VPs of model 
C. The lower panels show the VPs for a model with the same BH mass and mass-to-light ratio as 
model C, but for which the modified approach is used. For R = 0.0" and R = 4.83" the model 
VPs of the two different approaches are almost identical. However, for the larger radii, the VPs 
of model C clearly reveal large amounts of counter rotation. The VP of the upper-right panel has 
such a large counter-rotating component, that the best-fitting Gaussian no longer corresponds to 
the peak at positive velocities: it is very broad and centered around F ~ kms^^. This causes 
the large difference between the model predictions and the observed values at the outermost WHT 
point in Figure ^ Although the reconstructed V and a may differ strongly from the observed 
values, the hi and /12 values are very similar. Since these are the values that enter into the 
computation of x^, this discrepancy does not affect the statistics. 

It is interesting that the central velocity dispersions, as measured with the WHT, can be fit 
without the requirement of a central BH. This contradicts the conclusions reached from the Jeans 
modeling, which suggests, on the basis of the WHT measurements alone, that a BH of ~ 3 x 10* 
is required. The three-integral modeling, however, shows that only the HST measurements are of 
sufficient spatial resolution to discriminate strongly between models with and without a massive 
BH (cf. panels b and c of Figure 0). Clearly, without these high spatial-resolution kinematics, 
the case for a BH is only marginal: the plot in panel b of Figure |^ suggests a BH mass of 
~ 2.0 X 10* Mq, but cannot rule out models without a BH to a significant confidence level. 



7.2. Alternatives to a black hole 
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7.2.1. A dense cluster 

Although the dynamical evidence for the presence of a MDO of a few 10^ in NGC 4342 is 
compelling, it does not automatically imply evidence for a BH. Alternatives to a point mass, such 
as a cluster of brown dwarfs or stellar remnants are not ruled out by the modeling presented above. 
Any of these alternatives is only viable if its lifetime is not significantly smaller than the age of 
the galaxy (typically ~ 10 Gyr). There are three main processes that determine the evolution of 
a dense cluster: (i) core collapse, (ii) evaporation due to weak gravitational scattering, and (iii) 
physical collisions between the objects comprising the cluster. The latter of these processes is 
likely to ultimately lead to the formation of a single dense object, probably a massive BH. After 
a timescale TcoU, each object has physically collided with another (in a statistical sense). TcoW 
strongly depends on the mass and density of the cluster, as well as on the mass and size of the 
constituents (see Maoz 1997 and references therein). The characteristic timescale for evaporation 
of a cluster which consists of equal mass constituents is Tcvap ~ 300 Tj-dax (Spitzer & Thuan 1972), 
whereas typically after t^c ~ IBrreiax! a Plummer sphere of equal mass objects undergoes core 
collapse (Cohn 1980). Here Tj-eiax is the median relaxation time (see Spitzer & Hart 1971). By 
choosing a Plummer model for the dark cluster we are conservative, in that more concentrated 
clusters have shorter collapse times, and are thus less likely as alternatives for a massive BH. 

In order to constrain the size of a dark cluster, we construct models in which we replace 
the point-mass potential of the BH by a Plummer potential with a scale length e. We consider 
model B for which we replace the 3.6 x 10* Mq BH by a Plummer potential with the same mass. 



but with different values of e. The stellar mass-to-light ratio is kept constant at 6.25. Figure 10 
shows the resulting as function of e. The best fit is obtained for e = 0.0 (model B), and the fit 
deteriorates with increasing scalelength of the Plummer potential. The dotted line indicates the 
formal 99.73 per cent confidence level, and at this level of confidence we can rule out dark clusters 
with e > 0.07". This upper limit on the scale length of the Plummer sphere corresponds to 5.1 pc 
at the assumed distance of 15 Mpc, implying a central density of the cluster > 6.7 x 10^ Mqpc"^. 

We calculate the relaxation timescale Treiaxj and the collision timescale Tcou in which we adopt 



the mass-radius relation for non-collapsed objects used by Goodman & Lee (1989). In Figure |ll 
we plot the characteristic timescales for core-collapse and collisional destruction of our dark 
cluster, as function of the mass ra of the cluster's constituents. As can be seen, the timescales for 
core collapse {t^c ~ 7 x 10^^ (m/ M©)"^ yr) and evaporation (revap ~ 20rcc) for the inferred dark 
cluster in NGC 4342 exceed the Hubble time for m <^ few Mq, and neither of these processes 
thus allows us to rule out a dark cluster as an alternative to a BH. The collision timescale TcoU 
is the most restricting, and we are close to being able to rule out dark clusters of non-collapsed 
objects with masses less than ~ 0.001 Mq. However, for clusters of brown dwarfs with masses of 
~ 0.08 Mq, Tcoii is still of the order of 10^^ yr, and such clusters can clearly not be ruled out by 
current observations. Clusters of collapsed stellar remnants, such as white dwarfs or neutron stars, 
have collision timescales that are even much larger than the value of TcoH plotted in Figure 
This is due to the much smaller collisional cross-sections of these collapsed objects. In fact, the 
collapse time for a dark cluster can be made arbitrarily long by giving its objects an arbitrarily 
small mass, and the problems with collisions and mergers can be avoided by assuming the cluster 
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to be a collisioiilcss gas of elementary particles. However, another cffeet might allow one to rule 
out such a dark cluster: the trapping of stars by the cluster due to dynamical friction. If the time 
scale for this process is smaller than the age of the galaxy, enouh stars will get trapped such that 
the cluster can no longer be considered dark. Alternatively, these trapped stars can merge and 
form a massive BH, thus spoiling the goal for which the cluster was introduced. The time scale for 
the trapping to occur is independent of the mass m of the cluster constituents, but does depend 
on the mass m* of the stars being captured, and is given by {m/m*)Tcc (as long as m* ^ m, 
see Quinlan 1996). For NGC 4342 we thus find that this capture time scale is only smaller than 
~ 10-'^'^ yr for stars with m* ^ 70 M0. We can thus not use this trapping- mechanism to put an 
appreciable constraint on the nature of a possible dark cluster in NGC 4342. 

In conclusion, it is clear that we are still a long way from being able to rule out dark clusters 
as viable alternatives for a massive BH in the center of NGC 4342. 



7.2.2. An end-on bar 



It has been argued that a bar observed end-on may mimic the presence of a BH (Gerhard 
1988). In particular, the nuclear disk in NGC 4342 may in fact be a (very thin) nuclear bar. 
However, there are several reasons why this interpretation is unlikely. First of all, axisymmetric 
models without BH fail to fit the central velocity dispersion, as measured with the HST/FOS, 
by ~ 130 kms-^ It seems unlikely that the elongated orbits in a bar potential can fix this. 
Furthermore, there are no indications that (the center of) NGC 4342 is triaxial. Although van den 
Bosch &; Emsellem (1998) have provided evidence that NGC 4570, a galaxy with a double-disk 
morphology similar to NGC 4342, has been shaped under the influence of a rapidly tumbling 
bar-potential, none of the characteristics found in NGC 4570 that led to this conclusion are 
apparent in NGC 4342. In addition, no minor axis rotation is found, and the small amount of 
isophote twist observed (see BJM98) is limited to the outer region of the galaxy, and is more likely 
to be associated with a small warping of the outer disk, probably induced by the small companion 
at ~ 30" SE. Finally, the steep cusp of the bulge of NGC 4342 makes the bar hypothesis unlikely, 
since the pressure support from this strongly cusped bulge probably assures stability against bar 
formation (cf. van den Bosch &: de Zeeuw 1996). 



7.3. The dynamical structure of NGC 4342 

During the orbit integrations we store, in addition to the fractional time spent by orbit k in 
meridional cell n, t^, the time-weighted first and second order velocity moments averaged over cell 
n. With the solution for the orbit weights 7 we can compute the internal dynamical structure of 
the model averaged over each cell n. These are given by 

(^a)n = ^^^%^, (19) 
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and ^ ^ ^ 

Here we use the subscript a to indicate either of the cyhndrical coordinates R, (p or z. The 
cell-averaged velocity dispersions can be computed according to (cia)^ = (fa)„ — (fa)n- 

In order to suppress noise, we average (fa)n ^^'^ {'^a)n '^^^^ cells n that have the same radius 
r but different (r and being the standard spherical coordinates). Given the strongly flattened 
shape of NGC 4342 and its multi-component structure, we decided to split the galaxy in two 
parts: we determine the internal dynamical structure in two cones; one with half-opening angle of 
30 deg centered around the equatorial plane, including the nuclear and outer disks, and the other 
one with half-opening angle of 60 deg centered around the minor axis, representing the bulge of 
NGC 4342. 

We investigate the dynamical structure of models A, B, and C using the modified approach 
(see Section 6.1) and regularization as described in Section 5.5 in order to suppress noise. As 
a check, we compare the dynamical structures before and after regularization (with the amount 
of smoothing chosen such that Ax^ = 1). The regularized model is found to have a smoother 
dynamical structure, but the main features are similar for both cases indicating that our results 
are not too sensitive to the particular method and amount of regularization adopted. 

In Figure |l^ we plot the dynamics in the 'equatorial cone' as function of radius. The upper 
panels show the rms velocities for models A, B and C in the range between r = 0.1" and r = 12" 
(corresponding to the regime where we have kinematic constraints along the major axis). The 
middle panels display the ratio cTa/ctotai for the same radial interval (cr^otai ~ Sa^a)- Since there 
is no streaming motion in the radial and vertical directions, (<tr) and (cr^) are equal to their 
respective rms velocities. The difference between (f^) and cr^ reflects the streaming motion of the 
model. The lower panels show the ratios (Tz/o'r and (Ttp/an- The three models differ predominantly 
in the inner ~ 3.0", reflecting the differences in BH mass, but are very similar outside of 3.0": 
clearly (v^) dominates the dynamics outside this radius in accordance with the rapid rotation of 
the outer disk. The rapid change in dynamical structure (from azimuthally anisotropic to radially 
anisotropic) going from 3" to 12" reflects the strong increase of (V^/o")obs over this radial interval 
(see Figure ^: the dynamically cold outer disk, mainly built up of close-to-circular orbits with 
low o"^, becomes the dominant mass component. Over the same radial interval, the projected 
ellipticity increases from 0.4 to ~ 0.7, and this thus explains the observed correlation between 
ellipticity and lobs/Vmod of the best fitting isotropic Jeans model plotted in Figure At R ^ 8", 
<^(t>/(^R ~ 0.75, not too different from the value in the solar neighborhood where (J(f,/(j{i = 0.6 
(Dehnen 1998). 

Models B and C have cTz/fyR remarkably constant at ~ 0.9, and are thus not too different 
from two-integral models (for which this ratio is exactly 1.0). This is very different from the case 
without BH (model A) for which az/cTR is a strong function of radius R. This is the reason that 
the two-integral, isotropic Jeans model discussed in Section 4, could not fit the HST rotation 
velocities without BH, whereas model A can. 

The main change going from model A to model C, is a strong increase of cr^/o"^ in the 
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inner ~ 3". In these inner regions the circular velocities increase strongly with increasing BH 
mass. Nevertheless, all three models provide an almost equally good fit to the observed rotation 
velocities. We checked our solutions and found that going from model A to C, the ratio of 
— Lj;-orbits over +L^-orbits increases in the radial interval 0.5" < R < 3.0". This causes the net 
streaming motions of all three models to be roughly similar despite the large differences in circular 
velocities. Furthermore, it explains the strong increase in a^f, observed when going from model A 
to model C. 

The dynamics of the bulge (represented by the cone with half-opening angle of 60 deg centered 
around the minor axis) is presented in Figure We only plot the results out to -R = 3.0", 
corresponding to the radial interval where we have kinematic constraints along the minor axis. At 
R ^ 1" the models are again similar, being dominated by rotation, albeit to a much lesser extent 
than along the major axis. Model A is radially anisotropic in the central region (R ^ 1"), and 
aR decreases steadily with BH mass going from model A to model C. The radial anisotropy of 
the central region of the bulge of model A is required to fit the high central velocity dispersions 
observed as good as possible. Close to the equatorial plane, the steep rotation curve observed 
prevents model A from being too radially anisotropic (see also Section 7.4). Model B is again 
remarkably close to two-integral form, having OzI^r ~ 1.0. 



7.4. The influence of central radial anisotropy 

We now investigate to what extent radial anisotropy can influence the central kinematics of 
NGC 4342. We solve the orbit weights for the model with T/ = 7.25 M©/ L© and Mbh = using 
two different sets of constraints. The first one uses all constraints except the rotation velocities 
measured with the FOS. The second one excludes only the FOS velocity dispersions from the 



constraints. The results for both models are shown in Figure 14, where we plot the predictions for 
the rotation velocities and velocity dispersions of both models. We only plot the predictions for 
the FOS kinematics: both models have almost indistinguishable WHT kinematics. The lower two 
panels in Figure 14 show the dynamics of the resulting models; we plot both {v]^)^^'^ and {v'^)^^'^ 
as function of radius. These quantities are averaged between 6 = Odeg and 6 = 90 deg. 

The model with the FOS rotation velocities excluded from the constraints (solid lines) reveals 
the highest central velocity dispersion achievable without a BH, while still fitting the WHT 
measurements. The model predicts a central rotation curve that is too shallow to fit the observed 
rotation velocities, and is strongly radially anisotropic in the center. The model excluding the 
FOS velocity dispersions (dashed lines) shows a very steep central rotation curve. The model 
is strongly azimuthally anisotropic in the center, resulting in a much smaller central velocity 
dispersion. Clearly, a model without a central BH cannot simultaneously fit the rotation velocities 
and velocity dispersions measured with the FOS. 
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7.5. The density distribution in the outer region of NGC 4342 

The best-fitting mass-to- light ratio of Tj = 6.3Mq/Lq is unusually large for the stellar 
population of an early-type galaxy (e.g., van der Marel 1991). The mass-to-light ratio is mainly 
set by the large rotation velocities measured at the outside of NGC 4342 (along the major axis). 
The large value found may indicate that NGC 4342 is embedded in a massive dark halo. We 
did not attempt to include such a component for two reasons: First, we are mainly interested in 
trying to fit the central kinematics and to examine the mass of a possible BH. The characteristics 
of a possible dark halo will not affect the central dynamics, and therefore will not alter our main 
conclusion that NGC 4342 harbors a massive BH. It may however result in a different dynamical 
structure of the outer bulge and disk. Secondly, since we can fit all kinematics without having to 
infer a dark halo, including such a component in the models will merely lead to a large range of 
different models that can fit the observed kinematics equally well. Only kinematics at much larger 
radii can constrain the presence and characteristics of such a possible dark halo (cf. R97). 

Since mass-to-light ratio scales with distance d as d~^, another explanation for the large 
value of T/ derived might be an Tindcrcstimate of the distance. In order to bring the inferred 
mass-to-light ratio in accordance with the average value for early-type galaxies, NGC 4342 has 
to be at about twice the assumed distance of 15 Mpc. This is however hard to reconcile with 
the observed heliocentric velocity of NGC 4342 of only 714 km s"-*^ (de Vaucouleurs et al. 1976). 
Galaxies out to 15 Mpc from the center of the Virgo cluster can still experience a Virgocentric 
infall of a few hundred kms^"*^: our local group falls towards Virgo with 200 kms^-*^ (Tammann 
&: Sandage 1985). A distance of 30 Mpc for NGC 4342, however, would imply a Virgocentric 
infall velocity of ~ 1500 kms~^ (where we have taken Virgo to be at 1100 kms~^), which seems 
too large. So, the small recession velocity of NGC 4342 implies that it can not be too far behind 
Virgo such as to have an appreciable effect on the inferred mass-to-light ratio. 

Another problem related to the outer density distribution is that the mass model used (the 
MCE model) falls off as exp(— r^) at large radii. This is probably not correct, but since we only 
have photometry that is limited to a small radial extent, we cannot examine this in more detail. 
The main requirement is that we have a proper mass model for the major part of the galaxy. Our 
MGE model has a total luminosity of Lj = 3.57 x 10^ Lq. We have shown in Section 3, that for 
this luminosity we derive B — V = 1.09, in good agreement with the average value for early-type 
galaxies. This therefore suggests that our MGE model covers the major part of the galaxy, and 
that the galaxy does not extend far beyond the radial extent of our MGE model. In addition, 
as is the case with a dark halo, a change in the luminosity density profile in the outer parts of 
NGC 4342 will not affect our main conclusions regarding the presence of a central BH. 



7.6. Comparison with other BH detections 

There are now about a dozen BH candidates known. The best cases are the center of the 
MW, where a proper-motion study has revealed a 2.5 x 10^ M© BH (Eckart & Genzel 1997), and 
NGC 4258, where water masers are found to be in Keplerian rotation around a 3.6 x 10^ M.q BH 
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(Miyoshi ct al. 1995). Other strong evidence for the existence of massive BHs comes from the 
observation of the Fe Ka line at 6.4 kcV in a number of Scyfcrt 1 nuclei: the X-ray emission line 
exhibits relativistic motions, interpreted as arising from the accretion disk surrounding a massive 
BH (Tanaka et al. 1995; Nandra et al. 1997). 

Most other BH detections (or confirmations) are based on observations with the HST. 
Observations of nuclear gas disks with the HST has provided strong cases for massive BHs in 
M84, M87, NGC 4261, NGC 6251, and NGC 7052. All other cases are based on stellar dynamical 
evidence: M31, M32, NGC 3115, NGC 3377, NGC 3379, and NGC4486B. Of these only M 32, 
NGC 3379 and NGC 4342 discussed here have been confronted with three-integral modeling 
thus far. As discussed in Section 1, one can only hope to properly detect a BH if one can 
observe its kinematics inside tbh- All currently detected BHs have tbh > 0.1", and since all 
these galaxies, except the nucleus of the Milky Way, have now been observed with the HST, all 
current BH-candidates fulfill this requirement. Remarkable enough, most of these galaxies (M31, 
M32, M87, NGC 3115, and NGC 4594) were already considered BH candidates when the angular 
resolution of the observations was still an order of magnitude larger, and similar to the angular 
size of the radius of influence of the inferred BH (Rix 1993). The BH in NGC 4342 has a relatively 
small radius of influence, with an angular size of ~ 0.4". Indeed we have shown that only the 
HST/FOS observations are of sufficient spatial resolution to rule out models without a central BH. 

Kormendy & Richstone (1995) have suggested a relation between the mass of the BH and 
that of the bulge of its parent galaxy: (MBH/-^buige) = 0.22lQ;Qg x 10~^. For more recent 
discussions on this relation, including more recent BH detections see Kormendy et al. (1998), 
Ho (1998), van der Marel (1998), and van der Marel & van den Bosch (1998). NGC 4342 is the 
galaxy with currently the second-highest ratio of BH mass over bulge mass of 2.6^J| X 10-2, and 
is only superseded by the peculiar galaxy NGC 4486B, for which the case of a massive BH is 
less strong (Kormendy et al. 1997). NGC 3115 (MeH/Mbuige = 2.4 x 10"^) and the Milky Way 
(MenZ-^buige = 0.017 x 10"^) are other clear deviants. Clearly, the scatter around the suggested 
relation is considerable, seemingly as large as nearly two orders of magnitude. 



8. Conclusions 

Spectra obtained with the WHT and HST/FOS of the edge-on SO galaxy NGC 4342 have 
revealed a very steep central rotation curve and a strong central increase in velocity dispersion. 
These data suggest a large central mass concentration. In this paper we presented detailed 
dynamical models of NGC 4342 used to investigate whether its nucleus harbors a massive BH. 

We model the luminous density distribution of NGC 4342 with multiple Gaussian components. 
After projection and PSF convolution this model provides an excellent fit to the HST 7-band 
surface brightness distribution. The parameters of this model were derived with the MGE method. 

Simple isotropic Jeans models suggest that NGC 4342 harbors a massive BH of a few times 

10^ Mq. The actual mass of the BH depends on the data-set fitted: the WHT data suggest 
Mbh « 3 X 10^ Mq; the HST/FOS data suggest a somewhat larger BH mass of 6 x 10^ Mq. This 
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discrcpancy already suggests that the assumptions underlying the Jeans models, i.e., / = f{E, L^) 
and therefore (tr = az, are incorrect. This is also evident from the fact that the Jeans models 
cannot accurately fit the major-axis rms velocities measured with the WHT. These rms velocities 
are independent of the freedom in the anisotropy a^/aji allowed in the Jeans modeling. We 
find that for a mass-to-light ratio T/ = 6.2 Mq/ Lq and an isotropic velocity distribution 
{aR = dz = cr^) the Jeans model with Mbh = 3 x 10^ Mq provides the best fit to the observed 
WHT velocity dispersions along the major axis. However, the rotation velocities are not very well 
fitted and we find a correlation between Fobs/Knod and the local ellipticity of the projected surface 
brightness, such that the model underpredicts the rotation velocities in the highly flattened regions 
(dominated by the disk light) and overpredicts them in the less flattened region (dominated by 
the bulge light). This suggests that the different components in NGC 4342 have different velocity 
anisotropies. 

We thus constructed three-integral axisymmetric models of NGC 4342 in order to examine 
the mass of a possible BH and the dynamical structure of the different components. The modeling 
technique is an extension of Schwarzschild's orbit-superposition technique, and is based on finding 
the ensemble of orbits that best fits the observations. These models make no assumption about 
the dynamical structure and are fully general. This technique, developed by Rix et al. (1997) 
and Cretton et al. (1998), has previously been used to prove the existence of a massive BH 
of (3.4 ± 0.7) X 10^ Mq in the compact elliptical M32 (van der Marel et al. 1998). We have 
constructed a range of dynamical (Mbh, 'r7)-models of NGC 4342 to determine a central BH 
mass of 3.0t|;o ^ ^^^d an /-band mass-to-light ratio of 6.3^0 4 M0/ L0. The high spatial 

resolution of the HST/FOS data allow us to rule out models without a BH to a confidence level 
better than 99.73 per cent. With a similar confidence we can rule out models with a BH more 
massive than 7 x 10^ Mq. This upper limit on the BH mass is mainly due to the VP shape 
parameters and /i4. With the current data we can not rule out alternatives to a massive 
BH, such as a cluster of brown dwarfs or stellar remnants. Nevertheless, the QSO paradigm, 
together with the fact that the presence of massive BHs in galactic nuclei has unambiguously been 
demonstrated in a few galaxies were the inferred central densities are high enough to rule out 
dark clusters as alternatives (see discussion in Maoz 1997), make the interpretation of the inferred 
MDO in NGC 4342 in terms of a massive BH the most likely. 

We computed the intrinsic mean velocities and velocity dispersions of the three-integral 
models. The dynamical structures of the best fitting models vary strongly with radius, reflecting 
the multi-component structure of NGC 4342. Between 2" and 12" in the equatorial plane the best 
fitting models change from azimuthally anisotropic to radially anisotropic, while az/aji ^ 0.9. 
This explains the correlation between the projected ellipticity and the failure of the isotropic 
Jeans models to fit the observed rotation velocities along the major axis. The bulge in the best 
fitting model without BH is radially anisotropic. However, we have shown that even without 
the constraints of the measured HST/FOS rotation velocities, models without BH cannot fit the 
central HST/FOS velocity dispersion. The rotation velocities measured from the ground already 
constrain the amount of central radial anisotropy such that models without a BH cannot fit the 
high central velocity dispersion measured with the HST/FOS. 

The BH mass thus derived contributes a fraction of 2.6jlQ;| per cent to the total mass of the 
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bulge (1.2 X 10^° Mq). With this BH mass, NGC 4342 has one of the highest ratios of BH mass 
over bulge mass. Currently, the BH in our own galaxy has, with 0.02 per cent, the lowest BH 
mass to bulge mass ratio known: the scatter in the Mbh vs. Mbuige relation seems to be as large 
as two orders of magnitude. Extremely high spatial resolution is required in order to investigate if 
other galaxies have even lower values of MBH/-^buige- In conclusion, current data are consistent 
with a relation between bulge mass and BH mass, but the scatter is very large, and it is likely 
that the current MBn/^^bulgc ratios found are biased towards an upper limit. Although the newly 
installed Space Telescope Imaging Spectrograph (STIS) is likely to detect many more BH cases in 
the coming years, detection of BHs with masses of the order of a 0.02 per cent of the bulge mass or 
less in galaxies in Virgo or beyond, will probably have to await a next generation space telescope. 

Part of the calculations presented in this paper were performed on two UltraSparcs, kindly 
made availably by the Lorentz Center. Many thanks arc due to Eric Emsellem, for his help with the 
MGE analysis, and to Tim de Zeeuw, Walter Jaffe, HongSheng Zhao and Roeland van der Marel 
for many fruitful discussions. FvdB was supported by a Hubble Fellowship, #HF-01102.11-97A, 
awarded by STScI. NC acknowledges the hospitality of Steward Observatory where part of this 
work was done. 
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Fig. 1.— Contour maps of the WFPC2 /-band image of NGC 4342 at two different scales: 32" x 32" 
(left-hand panel) and 8" x 8" (right-hand panel). Superimposed are the contours of the MGE model 
of the intrinsic surface brightness convolved with the HST PSF (see Section 3). 

Fig. 2. — Observed rotation velocities V and velocity dispersions a (as determined from the best- 
fitting Gaussian, see text) along the major axis of NGC 4342 obtained with the WHT (crosses) 
and the FOS (solid dots). The gradient of the rotation velocity and the central velocity dispersion 
increase considerably going to the four times higher spatial resolution of the FOS. See BJM98 for 
details on the data. 

Fig. 3. — Results of the Jeans modeling. The solid dots with errorbars indicate the observed rotation 
velocities and velocity dispersions. Overplotted are four models that differ only in the mass of the 
central BH (0, 3, 5 and 10 x 10^ Mq). All models have fully isotropic velocity dispersions, and a 
stellar mass-to-light ratio of T/ = 6.2M0/LQ. The three panels on the left show the WHT kinematics 
and the model predictions for Vrot, <^p and Vrms- The model with Mbh = 3 x 10*^ M0 provides the 
best fit to the velocity dispersions. Neither of the four models provides an accurate fit to the 
rotation velocities (see also Figure^. The panels on the right compare the model predictions with 
the HST/FOS kinematics. Both the rotation velocities and the velocity dispersions suggest the 
presence of a BH with a mass of ~ 6 x 10*^ Mq. 

Fig. 4. — The ratio Vobs/^od of the observed rotation velocities over the rotation velocities 
predicted by the isotropic Jeans model with Tj = 6.2Mq/L0 and Mbh = 3x 10^ Mq, as function of 
the local, projected ellipticity of the isophotes. Only results beyond 2" are shown. At smaller radii 
seeing influences the observed quantities significantly. There is a clear correlation in that the Jeans 
model underpredicts the rotation velocities in the moderately flattened region, and overpredicts 
the rotation velocities in the radial interval where the isophotes are strongly flattened, and thus 
dominated by the light of the outer disk component. 

Fig. 5. — The Gauss-Hermite moments {mn = 0, 1, 2) of a Gaussian VP expanded around another 
Gaussian with the same dispersion a. Results are plotted as functions of the velocity difference /S.V 
between the two Gaussians, expressed in units of a. For Al/ = 0, the two Gaussians are identical 
and /iQ = 1 and hm = (m = 1, 2, 3, ...). 

Fig. 6. — The ratio V/a as function of radius of the major-axis kinematics of NGC 4342 as measured 
with the WHT (solid dots) and the HST/FOS (open circles). The dotted line indicates outside 
which V/a > 1.5. In the modified approach, we exclude counter-rotating orbits with circular radii 
Rc > Riim from the NNLS. 

Fig. 7. — Contour plots of = Xobs + ^sc which measures the goodness-of-fit to the constraints as 
function of BH mass, Mbh, and /-band mass-to-light ratio, T/. The first three contours define the 
formal 68.3, 95.4 and 99.73 per cent confidence levels (latter one is plotted with a thick contour). 
The subsequent contours are characterized by a factor two increase in Ax^. Solid dots indicate 
actual model calculations. The surface of panel a results from excluding /13 and as well as 
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all HST/FOS measurements from the constraints. In panel b we have only excluded the HST/FOS 
data, and panel c shows the resulting plot when all constraints are taken into account. The 
asterisks labeled A, B and C indicate special models discussed in the text. Finally, in panel d 
we have used the modified approach and included the additional constraint that counter-rotating 
orbits with Rc > 4.5" are not allowed in the solution. 

Fig. 8. — Kinematics for the three data sets (solid dots with errorbars) compared to the predictions 
for models A (dotted lines), B (solid lines) and C (dashed lines). Although we did not measure 
/13 and /14 from our HST/FOS spectra, we plot the model predictions in the lower two panels on 
the right. Model A, which has no BH, strongly underestimates the central velocity dispersion as 
measured with the FOS. The strong deviation of the predictions of model C for the outermost WHT- 
point along the major axis, is due to the 'artificial' counter-rotation of the model (see Figure 

Fig. 9. — Velocity profiles at four different radii along the major axis of NGC 4342. The observed 
VPs, reconstructed from the measured V, a, /13 and while assuming hm = for m > 5, are 
plotted as solid dots. The solid lines correspond to the model VPs, and the dashed lines to the 
Gauss-Hermite series (up to order 4) fitted to these model VPs. Upper panels correspond to model 
C, whereas lower panels correspond to a model with the same BH mass and mass-to- light ratio, but 
whose orbit solution was computed with the modified approach. The upper panels clearly reveal 
that the model VPs at larger radii are double peaked. This can yield strange Gaussian fits, as 
evident from the upper right panel. The modified approach excludes counter-rotating orbits with 
Rc > 4.5" from the orbit library, and does not reveal these double peaked VPs. 

Fig. 10. — The value of = Xobs + Xsc; which measures the goodness-of-fit to the constraints, as 
function of the scale length e (in arcsec) of the Plummer potential with a total mass of 3.6 X 10^ M0 
representing a cluster of dark objects in the center of NGC 4342. All models have a stellar mass- 
to-light ratio of 6.25. The best fit to the data is achieved in the limit e — > 0, which corresponds to 
our model with a BH rather than a dark cluster. We can rule out models with e > 0.07" at the 
99.73 per cent confidence level. 

Fig. 11. — The timescales for core collapse (tcc) and collisional destruction (tcou) of a 3.6 x 10^ 
Plummer cluster of non-collapsed objects with a central density of 6.7 x 10^ Mqpc~^ as function of 
the mass m of the constituents. If any of these timescales is less than ~ 10^*^ years, it makes such 
cluster an unlikely alternative to a BH. Unfortunately, current data does not allow us to rule out 
any of the dark clusters as an alternative MDO in NGC 4342. 

Fig. 12. — The dynamical structure of models A, B and C averaged over a cone with half-opening 
angle of 30deg centered on the equatorial plane. Upper panels show the rms velocities (w^)^''^ in 
kms~^ , middle panel the normalized velocity dispersions o"a/c7totah lower panels the ratios 
(^a/cTR- Solid curves are for the radial component (a = R), dotted curves for the azimuthal 
component {a = (j)), and dashed curves for the vertical component {a = z). Results are plotted over 
the radial interval where we have kinematic constraints along the major axis, i.e., 0.1" < R < 12". 



Fig. 13. — Same as Figure 12, except that here the quantities are averaged over a cone with 
half-opening angle of 60 deg centered on the symmetry axis i? = 0) of the models. 



Fig. 14. — Radial velocity anisotropy at work: Solid lines are for a model with T/ = 7.25 M©/ Lq 
and Mbh = 0, where we have neglected the FOS measurements of the rotation velocities. The 
dashed lines correspond to the same model, but for which we have neglected the FOS velocity- 
dispersion measurements. The two lower panels show the intrinsic dynamics, (v'^j)^^'^ (lower-left 
panel) and (vl)^''^ (lower-right panel), of the two models. 
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Table 1. Parameters of MGE model for the deconvolved /-band surface brightness. 



j 










(1) 


(2) 


(3) 


(4) 


(5) 


1 


490833.0 


0.032 


0.817 


1.40 X lO'^ 


2 


99417.9 


0.101 


0.865 


2.92 X 10'^ 


3 


67415.3 


0.282 


0.601 


1.07 X 10^ 


4 


84108.1 


0.343 


0.136 


4.47 X 10'^ 


5 


28511.8 


0.394 


0.856 


1.26 X 10^ 


6 


15529.3 


0.753 


0.622 


1.82 X 10^ 


7 


8490.8 


0.756 


1.000 


1.61 X 10^ 


8 


6055.4 


1.866 


0.665 


4.66 X 10^ 


9 


2951.4 


4.419 


0.250 


4.79 X 10^ 


10 


1572.8 


9.229 


0.266 


1.18 X 10^ 


11 


229.9 


11.854 


0.723 


7.76 X 10^ 



Note. — Column (1) gives the index number of each Gaussian. Column (2) gives its central 
surface brightness; column (3) its standard deviation (which expresses the size of the Gaussian 
along the major axis); column (4) its flattening; and column (5) its total 7-band luminosity. All 
Gaussians have the same position angle and the same center. 



